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Abstract 

We repeat our previous calculation of the spectrum distortion of massless neutrinos 
in the early universe with a considerably better accuracy and corrected for a missing 
numerical factor in one of the two ways of calculations presented in our paper 0] . Now 
both ways of calculations are in perfect agreement and we essentially reproduce our 
old results presented in the abstract of the paper and used in the calculations of light 
element abundances. We disagree with the criticism of our calculations presented in 
ref. [|. 
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1 Introduction 



In ref. Q we accurately calculated the distortion of the spectrum of primordial massless 
neutrinos due to the interaction with hotter electrons and positrons in the primeval cosmic 
plasma. Our results are in a reasonable agreement with the previous approximate or less 
accurate calculations (the list of references can be found in ref. [P). Recently a considerably 
different result was presented in ref. In view of that, we have redone our calculations 
with a better accuracy and also corrected for a missing numerical factor in one of the two 
different but equivalent ways of calculations presented in ref. ||T[ (see below). Due to this 
error the results obtained by the two different methods were somewhat different. Now both 
ways of calculations are in perfect agreement and we essentially reproduced our old results 
with a better precision]^. 

The increase in the energy density of neutrinos and the distortion of the spectrum of elec- 
tronic neutrinos have little effect on light element abundance, but it was recently pointed 
out that the corresponding increased fraction of cosmic relativistic matter has an impact 
on the CMB anisotropics, which may be detectable by coming satellite experiments ^. 

2 Description of calculation procedure 
2.1 Two ways of calculation 

In our previous paper on the subject Q we used two different approaches for calculating 
the distortion of the neutrino spectra. We numerically solved the system of equations both 
for the full neutrino distribution functions, f^., and for the deviations from equilibrium, 
~ {fvj ~ fv'^)/ ft^- III the last case the contributions to the collision integrals for all 
the processes vanish for vanishing 5, except for the interactions of neutrinos with electrons, 

where the "driving force" term, proportional to the temperature difference between v and 

^The results quoted in the abstract of ref. [^] as well as the results used in the calculation of AYhe are 
the correct ones. 
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1.399077 
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Table 1: Two ways of calculation. 

e^, gives a nonzero contribution. Moreover, for small x < 1 this contribution is also very 
small because the temperature difference is approximately T-^jTy ~ 1 -|- 0.005x^'^ which is 
close to 1 for small x. Here we have defined T^, = 1/a. 

However, in the case with the full distribution we erroneously omitted the factor 1.22 (from 
the Planck mass) for some of the reactions, which created somewhat different results for the 
two different approaches. Now, after we corrected for this factor, the results are in perfect 
agreement. 

In the recent paper |^ the factor 2 in the rate of the reaction VoV^ '^a'^a is missing. 
Fortunately, the difference in neutrino energy density due to this factor is very small, ~ 10~^, 
which is well below the numerical precision. 

We can define the effective number of neutrino species at asymptotically large time as: 

"eff = ^^4^^. (I) 

Pu p-y 

where the photon energy density is pj = (7r^/15)(aT^)^ and the equilibrium quantities are 
pf} = (7/8)(^Vl5)/«' and = {tt^ 15) {aT^i)\ 

In the last two columns in tables 1 and 2 we present N^q from eq. (|l]). The error on N^q, 
SNqQ, comes from the numerical error in the definition of aT^"^ = 1.40102 (see section p. 41) . 
We denote this error 6{aT^^) in tables 1 and 2. 
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2.2 Initial conditions 



Now, let us discuss the choice of the initial time Xj„. We made the runs for the system 
of kinetic equations with three different values Xj„ = 0.1,0.2 and 0.5. We found that the 
results of the runs with Xj„, = 0.1 and x^^ = 0.2 are the same with an accuracy of 10^^. 
This means that for x < 0.2 we can neglect the non-equilibrium corrections to the neutrino 
distribution functions. 

Let us note that already at Xin = 0.1 the dimensionless photon temperature differs from 
unity, aT ^ 1. For our calculations we took two possible sets of initial conditions. The first 
was used in the paper |||]: 

= feg = , aT{xin = 0.1) = 1.00006 . (2) 

These conditions correspond to separate energy conservation in the electromagnetic plasma 
before the time Note, that even though aT is very close to 1, we need to keep aT ^ 1, 
because our precision for the equilibrium temperature is of the order 0.00003. The second 
set of initial conditions is similar to the one we used in the paper Q: 



= ^y/T^^ ' ar(xi„ = 0.1) = 1.00003 . (3) 

These conditions correspond to the neutrinos being in thermal equilibrium with the electro- 
magnetic plasma before We found that both initial conditions give the same results for 
the neutrino energy density and other essential quantities for Xin = 0.1, i.e. the difference 
in the results is less than our numerical errors. If one instead would choose Xm = 0.2 then 
the condition (^) is more precise. 

2.3 Momentum grid 

We took the dimensionless momentum interval < y < 20. For equilibrium neutrinos 
p^{y > 2Q)/ py « 3 X 10^^. Because the nonequilibrium correction to the neutrino energy 
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density is of the order 1%, neglecting y > 20 does not affect the result for the neutrino 
energy density even if the calculation is done with 1% precision of the effect. There is a 
somewhat bigger correction in the reaction rates due to the preexponential factor, ~ p^, 
but even for that the momentum cut-off at y = 20 provides a sufficient accuracy. 

In order to choose the distribution of the momentum grid properly let us take a look on 
the differential energy density of the neutrinos dpy/dy = {1 / T['^)y^ f [y) . 97.5% of the energy 
density comes from particles with momentum in the interval 1 < y < 10. Particles with 
momentum 0.1 < y < 1 give 1.4% of the total energy density and particles with 10 < y < 
20 give 1.1% of the total energy density. Also note that non-equilibrium corrections are 
particularly important for particles with large momenta. 

These arguments advocate the use of a linear distribution in the region < y < 20 or 
log(y) distribution in the region 0.1 < y < 20. We found that the difference between these 
approaches with the same number of points in grid is about 10~^ for the neutrino energy 
density. 

The authors of the paper Q chose log(y) distribution of points in grid 10^^'^ < y < 10^'^ 
with 40 points per decade. With such a choice more than half of the points lie in the region 
y < 0.1, which gives only 0.0002% contribution to the neutrino energy density. In the most 
important decade 1 < y < 10 they have only 40 points. 

In order to check the errors connected with a finite number of points in grid we took the 
100, 200 and 400 point grids. The results hereof are presented in table 1. 

2.4 Time evolution 

We used three different methods of time evolution. 

1) Euler method. We control the errors connected with a finite number of points in time 
X in the following way. First, we run the program with some fixed number of points in x, 
distributed in the time interval Xj„ < x < xj in such a way that the distribution functions 
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do not change significantly at any momentum point y during one time iteration dx. Then 
we run the program for the entropy conservation law (i.e. with equilibrium neutrinos) with 
the same values of time Xj as in the first run. Finally we compare the asymptotical values 
of the temperature ratios with the theoretical value which is T^/Ty = (11/4)^/'^ = 1.40102. 
In order to have good precision we require that the numerical error in these temperature 
ratio should not be larger than ^ 0.00003 (fourth column in table 1). 

2) Bulirsch-Stoer method. Instead of the simple time evolution we used the Bulirsch-Stoer 
method, described in the book 

3) Method for stiff equations. In order to compare our results with the results of the 
paper we made calculations with their method for stiff equations. 

We found that the most precise (in calculation of T^'^) is the Bulirsch-Stoer method, but in 
the region of small time Xj„ < x < 1 it requires 10 times more processor time than the Euler 
or the stiff method. Moreover, in the case that we take Xj„ ^0.1 only the stiff method 
takes a reasonably small number of time steps. 

In the region of large time x > 1 the situation is the opposite. The Bulirsch-Stoer method 
requires 10 times fewer time steps than the Euler method. Unfortunately, we cannot control 
the precision of the stiff equations method in this time region. The problem is that the 
energy conservation law, which we use for the evaluation of the photon temperature, is not 
a stiff equation. In our calculations with the stiff method we therefore evolved the photon 
temperature as = T^'' + 6T. The value of the equilibrium photon temperature T^'^ is 
taken from the entropy conservation equation, while for the evolution of the small 6T we 
use the same time steps as we have for the kinetic equations in the stiff method. 

In the calculations with the corrections, 5{x,y), to the distribution functions we used the 
Euler method. In calculations with the total distribution functions, f{x,y), we used the 
Euler method for x < 1 and the Bulirsch-Stoer method for x > 1. In order to compare 
our results with the paper ^ , we also made the calculations using the stiff method in two 
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f{x,y) 
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Table 2: Different time evolution algorithms. 



different ways: with the stiff method evolution for all x, and with a combination of the 
stiff method for x < 1 and the Bulirsch-Stoer method for x > 1. In table 2 we compare 
the results for these ways of calculation with 100 points grid. We found that all ways of 
calculation give an effective number of neutrino species around N^ff = 3.0340, except the 
stiff method which gives a slightly larger value Ngff = 3.0344. 

3 Conclusion 

We have seen, that the effect of non-equilibrium neutrinos can be calculated with a very 
good accuracy if one takes a large enough number of momentum points in the important 
region 1 < y < 10 and a precise algorithm for the time evolution. Finally let us present our 
results for late times. Photon temperature: {aT^) final = 1.39910 it 0.00003 , correction to 
the energy density of electron neutrino: 5pi,^/ p^^ = 0.946 it 0.001 , correction to the energy 
density of muon neutrino: 8pi,^/ pi,^ = 0.398 it 0.001 , and effective number of neutrino 
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species: N^ff = 3.0340 ± 0.0003 
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